/*
 * Source code copied form test_s2_semi_memo_for.c
 * to do the FORWARD spherical harmonic transform.
 * output spherical harmonic descriptors of an obj file
 * run as: ./sphdes [#_out] [bw]
 */

#include <errno.h>
#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <vector>
#include <string>


#include "fftw3.h"
#include "makeweights.h"
#include "cospmls.h"
#include "FST_semi_memo.h"

using namespace std;

int main(int argc, char **argv)
{
  FILE *fp ;
  int i, bw, size ;
  int l, m, dummy;
  int cutoff, order ;
  int rank, howmany_rank ;
  double *rdata, *idata ;
  double *rcoeffs, *icoeffs ;
  double *weights;
  double *seminaive_naive_tablespace, *workspace;
  double **seminaive_naive_table ;
  fftw_plan dctPlan, fftPlan ;
  fftw_iodim dims[1], howmany_dims[1];

//  if (argc < 4)
//    {
//      cout<<"Usage: test_s2_semi_memo_for sampleFile outputFile bw [output_format]"<<endl;
//      exit(0);
//    }

  bw = atoi(argv[2]);

  /*** ASSUMING WILL SEMINAIVE ALL ORDERS ***/
  cutoff = bw ;
  size = 2*bw;

  /* allocate memory */
  rdata = (double *) malloc(sizeof(double) * (size * size));
  idata = (double *) malloc(sizeof(double) * (size * size));
  rcoeffs = (double *) malloc(sizeof(double) * (bw * bw));
  icoeffs = (double *) malloc(sizeof(double) * (bw * bw));
  weights = (double *) malloc(sizeof(double) * 4 * bw);
  seminaive_naive_tablespace =
    (double *) malloc(sizeof(double) *
		      (Reduced_Naive_TableSize(bw,cutoff) +
		       Reduced_SpharmonicTableSize(bw,cutoff)));
  workspace = (double *) malloc(sizeof(double) * 
				((8 * (bw*bw)) + 
				 (7 * bw)));


  /****
       At this point, check to see if all the memory has been
       allocated. If it has not, there's no point in going further.
  ****/

  if ( (rdata == NULL) || (idata == NULL) ||
       (rcoeffs == NULL) || (icoeffs == NULL) ||
       (seminaive_naive_tablespace == NULL) ||
       (workspace == NULL) )
    {
      perror("Error in allocating memory");
      exit( 1 ) ;
    }

  /* now precompute the Legendres */
  cout<<"Generating seminaive_naive tables..."<<endl;
  seminaive_naive_table = SemiNaive_Naive_Pml_Table(bw, cutoff,
						    seminaive_naive_tablespace,
						    workspace);

  /* construct fftw plans */

  /* make DCT plan -> note that I will be using the GURU
     interface to execute these plans within the routines*/

  /* forward DCT */
  dctPlan = fftw_plan_r2r_1d( 2*bw, weights, rdata,
			      FFTW_REDFT10, FFTW_ESTIMATE ) ;
      
  /*
    fftw "preamble" ;
    note that this plan places the output in a transposed array
  */
  rank = 1 ;
  dims[0].n = 2*bw ;
  dims[0].is = 1 ;
  dims[0].os = 2*bw ;
  howmany_rank = 1 ;
  howmany_dims[0].n = 2*bw ;
  howmany_dims[0].is = 2*bw ;
  howmany_dims[0].os = 1 ;
  
  /* forward fft */
  fftPlan = fftw_plan_guru_split_dft( rank, dims,
				      howmany_rank, howmany_dims,
				      rdata, idata,
				      workspace, workspace+(4*bw*bw),
				      FFTW_ESTIMATE );


  /* now make the weights */
  makeweights( bw, weights );
for (int dirid=1; dirid<atoi(argv[1])+1; dirid++){
	for (int rindex=1; rindex<33; rindex++) {
		char filename[30];
		sprintf(filename, "../data%d/%d.dat",dirid,rindex);
		/* now read in samples */
		fp = fopen(filename,"r");
		for(i = 0 ; i < size*size ; i++ ) {
			/* first the real part of the sample */
			fscanf(fp, "%lf", rdata + i );
			/* now the imaginary part */
			fscanf(fp, "%lf", idata + i );
		}
		fclose( fp ) ;

  
		/* now do the forward spherical transform */
		FST_semi_memo(rdata, idata,
		rcoeffs, icoeffs,
		bw,
		seminaive_naive_table,
		workspace,
		0,
		cutoff,
		&dctPlan,
		&fftPlan,
		weights );
    
  cout<<"about to write out coefficients"<<endl;

  /* now write out coefficients, but in what format ? */
  if ( argc == 4 )
    order = atoi(argv[3]);
  else
    order = 0 ;

  double temp=0; 
  
  char ofilename[30];
  sprintf(ofilename, "../out%d/%d.dat",dirid,rindex);
  fp = fopen( ofilename, "w" );
  if ( order == 0 )    /* code format */
   // for( i = 0 ; i < bw*bw ; i ++ )
   //   fprintf(fp, "%.15f\n%.15f\n", rcoeffs[i], icoeffs[i]);
  for ( l = 0 ; l < bw ; l++ ) {
	temp=0;
    for ( m = -l ; m < l + 1 ; m++ ) {
	  dummy = seanindex(m, l, bw);
	  temp+=(pow(rcoeffs[dummy],2)+pow( icoeffs[dummy],2));
	}
	temp=sqrt(temp);
    fprintf(fp, "%.15f\n", temp);
  }
  else                  /* human format */
    for ( l = 0 ; l < bw ; l++ )
      for ( m = -l ; m < l + 1 ; m++ )
	{
	  dummy = seanindex(m, l, bw);
	  fprintf(fp, "l = %d\t m = %d\t %.15f + %.15f I\n",
		  l, m, rcoeffs[dummy], icoeffs[dummy]);
	}
  fclose(fp);

  	  cout<<"finished writing coefficients"<<endl;
	}
}
  /* clean up */
  fftw_destroy_plan( fftPlan );
  fftw_destroy_plan( dctPlan );

  free(workspace);
  free(seminaive_naive_table);
  free(seminaive_naive_tablespace);
  free(weights);
  free(icoeffs);
  free(rcoeffs);
  free(idata);
  free(rdata);

  return 0 ;
}

